* Replicates all the results presented in the main text

*************************************************************************
* Load data
use "d.dta", clear

*************************************************************************
* Figure 4: Case selection
est clear

glm numerofirme i.bands c.pop2021 c.densita_x_km2 c.old_11 c.Tassodidisoccupazione c.perc_man_01 c.perc_agr_01 c.share_university_cens2011 c.digital_divide_rf c.mountain_perc i.COD_PROV if pop2021 < 25000 & resistent_region == 1, f(nbinomial) vce(cluster COD_PROV)
predict res if pop2021 < 25000 & resistent_region == 1, ansc

* produce sample mean local
sum w_sign, det

* case selection plot
twoway (scatter res w_sign if bands == 1 & regione != "Toscana" & pop2021 < 25000 & resistent_region == 1, mcolor(black%10) msymbol(o) legend(off) yline(0, lcolor(black%10) lpattern(dash) lwidth(medthick)) ytitle(Residuals) xtitle(Signature Rate) xline(`r(mean)', lwidth(medthick) lcolor(black%10) lpattern(dash)) xlabel(0 "0" .05 "0.05" .10 "0.10" .15 "0.15", format(%9.1f)) text(18.5 0.007 "μ", color(black%10) size(medsmall))) 

graph export "case_selection.png", replace width(3900)	

*************************************************************************
* Figure 5: Baseline + Matching Estimates

est clear

* Baseline Results
eststo Baseline: glm numerofirme bands c.pop2021 c.densita_x_km2 c.old_11 c.Tassodidisoccupazione c.perc_man_01 c.perc_agr_01 c.share_university_cens2011 c.digital_divide_rf c.mountain_perc i.COD_PROV if pop2021 < 25000 & resistent_province == 1 , cl(COD_PROV) f(nbinomial)

* Inverse Probability Weighting
logit bands c.pop2021 c.densita_x_km2 c.old_11 c.Tassodidisoccupazione c.perc_man_01 c.perc_agr_01 c.share_university_cens2011 c.digital_divide_rf c.mountain_perc i.COD_PROV if pop2021 < 25000 & resistent_province == 1, cl(COD_PROV)
predict pscore

sum pscore if pscore != ., det
sum pscore if bands == 0 & pscore != .
sum pscore if bands == 1 & pscore != .

gen support = 0 if pscore != .
replace support = 1 if pscore > .0047418

gen ips = .
replace ips = 1/pscore if bands == 1
replace ips = 1/(1-pscore) if bands == 0

* IPW (trim upper and lower 5% of propensity score)
eststo IPW: glm numerofirme bands c.pop2021 c.densita_x_km2 c.old_11 c.Tassodidisoccupazione c.perc_man_01 c.perc_agr_01 c.share_university_cens2011 c.digital_divide_rf c.mountain_perc i.COD_PROV if pop2021 < 25000 & resistent_province == 1 & pscore > .05 & pscore < 0.95 [pweight=ips], f(nbinomial) cl(COD_PROV)


* Mahalanobis Matching
psmatch2 bands if pop2021 < 25000 & resistent_province == 1, outcome(numerofirme) mahal(pop2021 densita_x_km2 old_11 Tassodidisoccupazione perc_man_01 perc_agr_01 share_university_cens2011 digital_divide_rf mountain_perc latitude longitude) neighbor(1)

rename _weight w_mahal

tab provincia bands if w_mahal !=.

eststo MM: glm numerofirme bands c.pop2021 c.densita_x_km2 c.old_11 c.Tassodidisoccupazione c.perc_man_01 c.perc_agr_01 c.share_university_cens2011 c.digital_divide_rf c.mountain_perc i.COD_PROV if w_mahal != . [pweight=w_mahal], f(nbinomial) cl(COD_PROV)

* Entropy balance
ebalance bands pop2021 densita_x_km2 old_11 Tassodidisoccupazione perc_man_01 perc_agr_01 share_university_cens2011 digital_divide_rf mountain_perc latitude longitude if pop2021 < 25000 & resistent_province == 1

rename _webal w_eb

eststo EB: glm numerofirme bands c.pop2021 c.densita_x_km2 c.old_11 c.Tassodidisoccupazione c.perc_man_01 c.perc_agr_01 c.share_university_cens2011 c.digital_divide_rf c.mountain_perc i.COD_PROV if w_eb != . [pweight=w_eb], f(nbinomial) cl(COD_PROV)

coefplot (Baseline, label(Baseline Model)) (IPW, label(Inverse Probability Weighting)) (MM, label(Mahalanobis Matching)) (EB, label(Entropy Balancing)), keep(bands) xscale(range(0 0.5)) xlabel(0.0 0.1 0.2 0.3 0.4 0.5, format(%9.1f)) xtitle("Negative Binomial Coefficient") xline(0, lcolor(gs14) lpattern(dash) lwidth(medthick)) mlabel(@plot) mlabposition(1) legend(off)

graph export "matching.png", replace	

*************************************************************************
* Figure 6: Baseline + Pre-war Controls
est clear

eststo a: nbreg numerofirme bands c.pop2021 c.densita_x_km2 c.old_11 c.Tassodidisoccupazione c.perc_man_01 c.perc_agr_01 c.share_university_cens2011 c.digital_divide_rf c.mountain_perc i.COD_PROV if pop2021 < 25000 & resistent_region == 1, cl(COD_PROV) 

eststo b: nbreg numerofirme bands c.pop2021 c.psu1919_vv dlab1921 shcrop1921 ind_workers1921 literacy1911 sh_pop_1911_be6 c.mountain_perc i.COD_PROV if pop2021 < 25000 & resistent_region == 1, cl(COD_PROV)

coefplot (a, label(Baseline Model)) (b, label(Pre-war Controls)), keep(bands) xscale(range(0 0.5)) msymbol(O) xlabel(0.0 0.1 0.2 0.3 0.4 0.5, format(%9.1f)) xtitle("Negative Binomial Coefficient") xline(0, lcolor(gs14) lpattern(dash) lwidth(medthick)) mlabel(@plot) mlabposition(1) legend(off)

graph export "main.png", replace width(3900)

*************************************************************************
* Figure 7: Estimates by Political Leaning of the Local Band

gen PCIsoc = 0
replace PCIsoc = 1 if bands_d1_PCIsoc == 2

gen other_PCIsoc = 0
replace other_PCIsoc = 1 if bands_d1_PCIsoc == 1

label var PCIsoc "Communists & Socialists"
label var other_PCIsoc "Other Bands"

est clear
foreach bands in "PCIsoc other_PCIsoc" {
	
qui eststo: nbreg numerofirme `bands' c.pop2021 c.densita_x_km2 c.old_11 c.Tassodidisoccupazione c.perc_man_01 c.perc_agr_01 c.share_university_cens2011 c.digital_divide_rf c.mountain_perc i.COD_PROV if pop2021 < 25000 & resistent_region == 1, cl(COD_PROV)

}

coefplot est* , keep(*PCI*) labels msymbol(O) xscale(range(0 0.5)) xlabel(0.0 0.1 0.2 0.3 0.4 0.5 0.6, format(%9.1f)) legend(off) xtitle("Negative Binomial Coefficient")  xline(0, lcolor(gs14) lpattern(dash) lwidth(medthick)) coeflabels( , wrap(12)) 

graph export "ideology.png", replace	

*************************************************************************
* Figure 8a: Control for violence

est clear

eststo: nbreg numerofirme bands strage c.pop2021 c.densita_x_km2 c.old_11 c.Tassodidisoccupazione c.perc_man_01 c.perc_agr_01 c.share_university_cens2011 c.digital_divide_rf c.mountain_perc i.COD_PROV if pop2021 < 25000 & resistent_region == 1, cl(COD_PROV)

lincom bands - strage

coefplot est* , keep(bands strage) labels msymbol(O) xscale(range(0 0.5)) xlabel(0.0 0.1 0.2 0.3 0.4 0.5, format(%9.1f) labsize(large)) ylabel(,labsize(vlarge)) legend(off) xtitle("Negative Binomial Coefficient" , size(vlarge))  xline(0, lcolor(gs14) lpattern(dash) lwidth(medthick))

graph export "bandsVSviolence.png", replace	

*************************************************************************
* Figure 8b: Effect of violence by type of victim

est clear

eststo: nbreg numerofirme vittimetot c.pop2021 c.densita_x_km2 c.old_11 c.Tassodidisoccupazione c.perc_man_01 c.perc_agr_01 c.share_university_cens2011 c.digital_divide_rf c.mountain_perc i.COD_PROV if pop2021 < 25000 & resistent_region == 1, cl(COD_PROV)

eststo: nbreg numerofirme numero_vittime_civili numero_vittime_partigiani c.pop2021 c.densita_x_km2 c.old_11 c.Tassodidisoccupazione c.perc_man_01 c.perc_agr_01 c.share_university_cens2011 c.digital_divide_rf c.mountain_perc i.COD_PROV if pop2021 < 25000 & resistent_region == 1, cl(COD_PROV)

coefplot est* , keep(vittimetot numero_vittime_civili numero_vittime_partigiani) xscale(range(0 0.05)) ylabel(,labsize(vlarge)) coeflabels(vittimetot = `""All" "Casualties""'  numero_vittime_civili = `""Civilian" "Casualties""'   numero_vittime_partigiani = `""Partisan" "Casualties""' ) xlabel(0.00 0.01 0.02 0.03 0.04 0.05, format(%9.2f) labsize(large))  legend(off) xtitle("Negative Binomial Coefficient", size(vlarge))  xline(0, lcolor(gs14) lpattern(dash) lwidth(medthick))

graph export "typeofcasualties.png", replace	
